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ABSTRACT 

The migration of low mass planets, or type I planetary migration, has been studied 
in hydrodynamical disc models for more than three decades. For a long time, it was 
thought to be very rapid and directed inwards due to Lindblad torques. More recently, 
it has been shown that the corotation torque, linked to the horseshoe motion of the gas 
near the planet, may slow down or even reverse migration. How is this picture modified 
by the expected presence of a magnetic field in the protoplanetary disc ? When the 
magnetic field is strong enough to prevent horseshoe motion, the corotation torque is 
replaced by a torque arising from magnetic resonances which may significantly alter 
the migration rate. In the case of a weaker magnetic field, the magnetic field is not 
strong enough to prevent horseshoe motion and a corotation torque then exists. In this 
regime, recent turbulent MHD simulations have reported the existence of an additional 
component of the corotation torque due to the presence of the magnetic field. The aim 
of this paper is to investigate the physical origin and the properties of this additional 
corotation torque. 

We performed MHD simulations of a low mass planet embedded in a 2D laminar 
disc threaded by a weak toroidal magnetic field, where the effects of turbulence are 
modelled by a viscosity and a resistivity. We confirm that the interaction between the 
magnetic field and the horseshoe motion of the gas results in an additional corotation 
torque on the planet, which we dub the MHD torque excess. We demonstrate that it is 
caused by the accumulation of the magnetic field along the downstream separatrices 
of the horseshoe region, which gives rise to an azimuthally asymmetric underdense 
region at that location. The properties of the MHD torque excess are characterised 
by varying the slope of the density, temperature and magnetic field profiles, as well as 
the diffusion coefficients and the strength of the magnetic field. The sign of the torque 
excess and its radial distribution are found to be in agreement with the earlier full 
MRI simulations. This sign depends on the density and temperature gradients only 
and is positive for profiles expected in protoplanetary discs. The magnitude of the 
torque excess is in turn mainly determined by the strength of the magnetic field and 
the turbulent resistivity. It can be strong enough to reverse migration even when the 
magnetic pressure is less than one percent of the thermal pressure. The MHD torque 
excess can therefore lead to outward planetary migration in the radiatively efficient 
outer parts of protoplanetary discs, where the hydrodynamical corotation torque is 
too weak to prevent fast inward migration. 

Key words: planet-disc interaction - protoplanetary discs - accretion, accretion 
discs - magnetic fields - MHD 



1 INTRODUCTION 

The gravitational interaction between planets and their par- 
ent protoplanetary disc plays a prominent role in shaping 
planetary systems from the early stages of their formation 
and evolution. The torque exerted by the disc on a planet 



drives orbital migration, a change in the planets semi-major 
axis during the disc lifetime. The direction and speed of 
migration are intimately related to the planet mass, and 
the disc properties near the planet's orbital radius (tur- 
bulent viscosity, density and temperature profiles). Several 
regimes of migration are commonly distinguished, based on 
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the planet's ability to open an annular gap around its orbit. 
This ability depends on the planet to primary mass ratio, 
the aspect ratio of the disc h, and the viscosity parameter 
modelling the effect of turbulence, a. For typical disc aspect 
ratios and alpha viscosities in regions of planet formation 
(h ~ 0.05, a ~ a few xlO -3 ), type I migration applies to 
planets up to 10 to 20 Earth masses. Giant planets in the 
Jupiter-mass range carve a d eep gap around their orb it and 
experience type II migration ijLin fc Papaloizoulll986h . Sub- 
giant planets that open a partial gap may undergo type HI 
runaw ay migration in massive discs |Masset fc Papaloizou! 
2003). In this paper, we focus on type I-migrating planets. 

The torque acting on type I-migrating planets has two 
components: the differential Lindblad torque and the coro- 
tation torque. The differential Lindblad torque is the rate of 
angular momentum carried away by the spiral density waves 
(wakes) that the planet generates in the disc. Its sign and 
magnitude result from a sligh t asymmetry in the wakes' den- 
sity distribution (Ward 1997). For typical density and tem- 
perature profiles in discs, the differential Lindblad torque is a 
negative quantity. Alone, it would make the direction of type 
I migration inwards, on a time scale much sh orter than the 
typic al lifetime of protoplanetary discs (e.g., iTanaka et al.l 
l2002h . 

The corotation torque, or horseshoe drag, accounts 
for the angular momentum exchanged between the planet 
and the disc material within its coorbital region. It 
has been subject to intense investigation after it was 
discovered that its amplitude may be comparable to, 
or even exceed that of the differential Lindblad torque 
in radiatively inefficient regi o ns of planet formation 
(IPaardekooper fc Mellemal |200d: IBaruteau fc Massetl 20081 ; 



Paardckooper & Papaloizou 2008; Masset & Casoli 200 



— opt- _ 

Kiev et aljhoogi: IPaardekooper et al.ll2010l ; iMasset fc Casolil 
20ld ; iPaardekooper et al.l 2011 V The corotation torque is 



directly related to the distribution of the fluid's vortensity 
within the planet's horseshoe region (vortensity refers to the 
ratio of the vertical component of the vorticity to the surface 
density, and is also known as potential vorticity). In viscous 
disc models, the corotation torque is in general the sum of 
several terms: 

(i) A term proportional to the vortensity gradient across 
the horseshoe region, which arises from advection-diffusion 
of vortensity inside this re gion. It is o f ten called vor tensity- 
related corotation torque jwardlll99ll ; lMassedl200lh . 

(ii) A term proportional to the entropy gradient across 
the horseshoe region, similarly named entropy-related coro- 
tation torque. It originates from the (formally) singular pro- 
duction of vortensity at the downstream separatrices of the 
horseshoe region, location at which the entropy gets discon- 
tinuous due to the advect i on of e ntropy within the horsesho e 
region jMasset fc Casolil 120091 : IPaardekooper et alj I2010T) . 
In the limit of locally isothermal discs, where the radial 
profile of temperature is constant in time, the corota- 
tion torque features instead the local temperatur e gradient 
jCasoli fc Massed 120091 : IPaardekooper et alj|20ich . 

The sign and amplitude of the total corotation torque 
depend on the disc's density and temperature profiles across 
the planet's horseshoe region, as well as on viscosity and 
thermal diffusivity in this region. This is in contrast to the 
differential Lindblad torque, whose amplitude essentially de- 



pends on the local temperature gradient. The potentially 
large amplitude of the corotation torque makes it play a 
very important role in the evolution of young planetary sys- 
tems. In this context, much effort has been put forward to 
produce accura te, yet simple analytic formulae for the coro- 
tatio n torque (|Masset fc Casolil l2010l : IPaardekooper et "ail 
2011). They may help address some of the shortcomings 
recently spotted by models of planet population syntheses 
dlda fc Linl |200SI: iMordasini et alj 120091 ; ISchlaufman et all 
l2009l : lHellarv fc Nelsonll2012r i. 

For type-I migrating planets, the width of the horse- 
shoe region is a fraction of the disc's pressure scale height, 
and thus of the typical size of turbulent eddies. It is a priori 
unclear whether turbulence would generally act as a vis- 
cous and thermal diffusion over the horseshoe region of low- 
mass planets. Put another way, it is unclear whether the 
corotation torque in turbulent discs behaves similarly as in 
viscous discs. A few recent studies have started to exam- 
ine this issue |f3aruteau fc Linl 2010l: IBaruteau et alj l201ll : 
it al.ll2012T). Adopting the turbu - 



lUribe et al.ll201ll :l Pierens et al 



201 



lence model originally developed bv lLaughlin et all l|2004 ). 
where turbulence is generated by stochastic forc i ng, th e 2D 
hydrodynami c al sim ulations of IBaruteau fc Linl (|2010T l and 
IPierens ~et al. (2012) showed that both the differential Lind- 
blad torque and the corotation torque behave similarly as in 
equivalent laminar viscous discs. 

How is the picture modified by the expected presence of 
a magnetic field ? In the case of a strong azimuthal magnetic 
field that prevents horseshoe motion, the corotation torque 
is replaced by angular momentum carried away at mag- 
netic resonances, where the rotational velocity relative to the 
planet is equal to the propagation speed of a slow MHD wav e 
or an Alfven wave (|Terqueml |2003l ; iFromang et"afl 120051 ). 
In the intermediate case of a weak magnetic field, horse- 
shoe motion can be expected to take place and give rise 
to a corotation torque. This has been confirmed by simu- 
lations of discs in which the magneto-rotational instability 
(MRI) drives magneto -hydrodynamical (MHD) turbulence. 
IBaruteau et alj (|201lh (hereafter BFNM11) carried out 3D 
MHD simulations of discs with turbulence due to the MRI 
operating throughout the whole disc. They adopted a disc 
model with a weak mean toroidal magnetic field, in which 
non-ideal MHD effects and vertical stratification were ne- 
glected. They showed the existence of horseshoe dynamics 
and an unsaturated corotation torque for Saturn-like planets 
embedded in thick (h — 0.1) discs. They also found the exis- 
tence of an additional coro tation torque with moderate am- 
plitude. lUribe et alj (|2011 ) performed 3D MHD simulations 



Bar uteau et alj (|20 1 lT ) . except for 



with a similar setup as in 
the inclusion of vertical stratification. Their simulations with 
intermediate-mass planets showed outward type I migration, 
which is indicative of the presence of a fairly large and un- 
saturated positive corotation torque in their simulations. 

The aim of this paper is to investigate the physical ori- 
gin and the properties of this additional corotation torque 
due to a weak magnetic field, discovered by MHD turbulent 
simulations. For this purpose, we adopt a two-dimensional 
laminar disc model threaded by a weak toroidal magnetic 
field, and in which turbulence is modelled by viscous and 
magnetic diffusion. This simplified setup is chosen for several 
purposes. First, the lower computational cost of performing 
2D simulations allows a wide exploration of the parameter 
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space, so far inaccessible to 3D turbulent simulations. Sec- 
ond, the simplicity of the model allows a deeper physical 
understanding, which should prove useful to interpret more 
complex turbulent simulations. This approach is not new but 
has been carried out in several other contexts (see Section [2] 
below). Modelling the effects of turbulence by diffusion co- 
efficients is obviously a strong assumption, which should be 
checked for the problem considered in this paper by com- 
paring with numerical simulations of turbulent discs (a first 
step in this direction is done in Section 5 of this paper). We 
note, however, that such an approach has been successful 
in studies of migration in hydrodynamical disc models as 
mentioned earlier in this introduction. 

Our simulation results confirm the existence of an addi- 
tional corotation torque acting on low-mass planets embed- 
ded in discs with a toroidal magnetic field. A wide explo- 
ration of the parameter space allows us to determine numer- 
ically how the sign and magnitude of this additional corota- 
tion torque depend on the disc properties (radial profiles of 
density, temperature and toroidal magnetic field, viscosity, 
magnetic diffusivity) . For typical density and temperature 
profiles, the additional corotation torque is positive. Its am- 
plitude has a very steep dependence on the magnetic field 
strength and diffusivity, and in some cases it may largely 
exceed the amplitude of the differential Lindblad torque. 

This paper is organised as follows. The physical model 
and numerical setup are presented in Section [2] The topol- 
ogy of the magnetic field within the horseshoe region is de- 
scribed in Section^ The properties (sign, magnitude) of the 
new corotation torque in weakly magnetised discs are exam- 
ined in Section U A detailed comparison with the results 
of the 3D MHD turbulent simulations of BFNMII follows 
in Section [5] Concluding remarks and future directions are 
given in Section [5] 



2 PHYSICAL MODEL AND NUMERICAL 
SETUP 

2.1 Disc model 

We consider a two-dimensional disc model, assuming a lo- 
cally isothermal equation of state (i.e. the temperature de- 
pends on radius but not on time). The gas is assumed to 
be non-selfgravitating. The disc is laminar, and the effects 
of MHD turbulence are modelled by effective diffusion coef- 
ficients: a viscosity and a resistivity. Such an approach has 
been used in th e past in a variety of cont exts. Since the 
seminal work of IShakura k, Sunvaevl |l973), hydrodynam- 
ical disc models with an effective turbulent viscosity have 
been widely used to study many different aspects of disc 
dynamics, including the problem of planet migration. An 
effective resistivity has been used for the st udy of mag- 
netic flux transport in acc retion discs (e.g. iLubow et al.l 
1994 iGuilet fc Ogilyidl2012ft and the launch ing of jets (e.g. 
Ferreira fc Pelletierlll995l ; IZanni et all 120071 '). The viscosity 



is parameterised by the usual a parameter : v = ac^/flK, 
where c s is the sound speed and £Ik is the keplerian an- 
gular velocity. The value of the magnetic diffusivity r\ is 
then determined by the assumed magnetic Prandtl num- 
ber defined as V = v/f), giving : r\ — a/Vc^/Q,K . The ef- 
fective resistivity of MHD turbulence driven by the MRI 



has bee n measured in numerica l simulations of local disc 
models jGuan fc Gammid 120091 ; iLesur fc Longarettil |2009| ; 
iFromang fc Stond 12009 ) . They showed that the magnetic 
Prandtl number of MHD turbulence is of order unity, though 
it could differ from 1 by a factor of a few depending on the 
magnetic field configuration. In most of this article, we there- 
fore assume V = 1, and relax this assumption in Sections l4.3 
and [5] As for the radial dependence of the diffusion coeffi- 
cients, either uniform alpha parameter or uniform diffusion 
coefficients (i.e. v and r\) are considered, giving essentially 
the same results if the value at the planet's orbital radius is 
the same. 

The basic equations solved are : 



dv 



0. 



p ( + v ■ Vv ) = — pV$ — Vp 



H (V x B) x B + V • T, 

Po 



3B 

-— = V x (v x B - 77V x B) 



pv 



Vv + (Vv) 



~(V-v)I 



(1) 



(2) 



(3) 



(4) 



where p is the density, v the velocity, $ the gravitational 
potential, p the pressure, B the magnetic field, T the viscous 
stress tensor, r\ the magnetic diffusivity, v the kinematic 
viscosity and I the unit tensor of second rank. 

A cylindrical polar coordinate system (r, ip, z) is used, 
the disc occupying the plane z — 0. We assume a simple 
initial magnetic field geometry, with only a toroidal compo- 
nent B v . This is motivated by the fact that in sheared MHD 
turbulence the azimuthal component of the magnetic field 
tends to be the strongest one. Furthermore, we initialise the 
calculation with power law profiles for the surface density E, 
temperature T and azimuthal magnetic field B v . The power 
law indices are defined as : 

d log B v 



dlogr 

dlogE 
d log r ' 

rilogT 
dlogr 



(•») 



(6) 



(7) 



The strength of the magnetic field is measured in terms 
of the plasma parameter /3, which is defined as the ratio of 
the thermal to the magnetic pressure: 

P = -Pth/iW = 2c 2 Jv 2 A , (8) 

where va = B/y/pop is the Alfven velocity, and po is the 
vacuum permeability. 

The azimuthal velocity is initialised with a profile cor- 
responding to a balance between the gravity of the star, the 
centrifugal force, the pressure gradient and the magnetic 
force. In can be written (in a non rotating frame): 

2(b + l)~ 



1+ (p+q+ 



P 



0) 
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where flu = \/ GM, / r 3 is the Keplerian angular frequency, 
h = H/r is the aspect ratio of the disc with H = c s /Qk the 
disc scaleheight. 

2.2 Planet 

A planet is introduced in the disc at the beginning of the 
simulations and is held in a fixed circular orbit, at a radius 
r p and azimuth tp p — 0. We work in the frame rotating with 
the planet at angular frequency fi p and therefore include the 
required indirect terms in the equation of motion. The planet 
potential is smoothed over a softening length e = O.QH(r p ) 
(though this fiducial value is varied in Section I4.4f) and is 
expressed as: 



GM V 



?>| 2 + € 



2 "I V2 ' 



(10) 



The torque exerted by the disc on the planet is com- 
puted with the classical expression : 



S ~ p rdrdip 
dip 



(11) 



Unless otherwise specified, the torque calculations presented 
in this paper include the planet's Hill sphere, whose radius 
is: 



M p 
3M, 



1/3 



(12) 



The effect of excluding the Hill sphere is studied in Sec- 
tion [4T4] The Hill sphere is also excluded in Section [5] for the 
comparison with BFNM11. Finally, note that all the torques 
and torque density distributions presented in this paper are 
per unit mass of the planet, although we then refer to it as 
"torque" rather than "specific torque" for conciseness. 

2.3 Normalisation, fiducial parameters and 
migration regime 

Without loss of generality, we normalise the radius, surface 
density and time such that r p = 1, S(r p ) = 1, and Q, p = 1. 
The magnetic field is then normalised such that /Uo = 1. 
The free parameters of the problem are then the planet-to- 
primary mass ratio M p /M*, the value at the planet orbital 
radius of the aspect ratio of the disc h, the plasma param- 
eter P, and the diffusion coefficients (determined by a and 
V), and the slopes indices of the magnetic field b, the sur- 
face density p, and temperature q. We use as fiducial pa- 
rameters the following values: M p = 2 x 10 -5 M* (around 
7 earth mass's for a Sun- like star), h — 0.05, /3 = 100, 
q = 5.10~ 3 , V = 1, b = — 1. For p and q, we consider 
two fiducial disc models (the same ones as in BFNM11) : 
model 1 has p — — 1/2 and q = — 1 (corresponding to a 
uniform aspect ratio of the disc). Model 2 has p — —3/2 
and q — 0, which corresponds to the special case where the 
corotation torque vanishes in a non magnetised disc, because 
the vortensity and the temperature are uniform. Most of the 
parameters are then varied around this fiducial set of values. 

Migration regime. The fiducial value of the planet mass 
falls into the regime of type I migration, since M p /(M*/i 3 ) = 
0.16. The relative strength of the magnetic field and horse- 
shoe motion may be measured with the ratio of the Alfven 



speed to the shear velocity at the separatrix of the horse- 
shoe region. With the half- width of the horseshoe region in a 
hydrodynamical case estimated as: x s ~ l.lr p w M p /(M*h), 
one obtains: 



"A 



,(r p - 



0.86 



M p /3 ' 



(13) 



The fiducial parameters give: va/v^{x s ) ~ 0.21. Since the 
Alfven speed is significantly smaller than the shear veloc- 
ity at the separatrix of the horseshoe region, the magnetic 
field is not expected to prevent horseshoe motion of the gas. 
The numerical results indeed show the presence of horseshoe 
dynamics. 

Magnetic res onances. In th e pre sence of an azimutha l 
magnetic field, iTerqueml (|2003T ) and iFromang et al.l (|2005t) 
showed that magnetic resonances can arise at the locations 
where the shear velocity equals the propagation velocity of 
slow MHD waves (one resonance on each side of the planet). 
If the shear is the unperturbed Keplerian shear, then the 
separation between the two magnetic resonances (i.e. the 
width of the region where the flow is sub-slow MHD) can be 
expressed as: 

AH 

Ar mr = — . (14) 

3^1 + /V2 

With our fiducial parameters, this gives Ar rm = 9.3 x 10~ 3 . 
In the highest resolution runs of the convergence study pre- 
sented in the Appendix, this region is resolved by 6 — 7 cells. 
Yet, we do not see any evidence for magnetic resonances. The 
very good convergence of the results shows that numerical 
diffusion plays a minor role compared to physical diffusion, 
and suggests that the magnetic resonances are absent for a 
physical reason rather than a numerical reason. 

Two physical processes may be the cause for the absence 
of magnetic resonances in our setup: applied resistive diffu- 
sion and horseshoe dynamics. Since the postulated magnetic 
resonances are close to the corotation radius for our fiducial 
magnetic field strength, diffusion may play an important 
role. The time to diffuse from one magnetic resonance to 
the other is during which a fluid element on a mag- 

netic resonance moves by A(j> = 16/9V /a(l + f3/2)~ 3 / 2 h ~ 
0.05 ~ h for our fiducial parameters. Given this small dis- 
placement during a diffusion time, diffusion is expected to 
have a significant impact on the dynamics at the magnetic 
resonances. 

If they were present, the postulated magnetic reso- 
nances would lie within the horseshoe region. The shear 
is, however, radically modified inside the horseshoe region, 
thus changing the conditions for the magnetic resonances 
to appear. The shear velocity can be expected to equal the 
slow MHD speed in a small region close to the stagnation 
point, however it is unclear whether slow MHD waves can be 
launched in these conditions, which are very different from 
those in the absence of horseshoe motion. In conclusion, be- 
cause of these two physical reasons we may expect horseshoe 
motion slightly modified by the presence of the magnetic 
field rather than magnetic resonances, which is confirmed 
by the numerical results. 

Timescales. The properties of the corotation torque are 
sensitive to the ratio of the diffusion timescale to the u-turn 
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or libration timescales. The diffusion timescale across the 
horseshoe region is: 



Tdiff = — ~ 1.21 



(15) 



which gives rata — 6.2 orbits at the planet's orbital radius 
for the fiducial parameters. The libration timescale is: 



87rr D 



Tlib 



8tt \M,h 



(16) 



giving rub — 61 orbits for the fiducial parameters. And fi- 
nally, the U-turn timescale is approximately ru-tum — hrub, 
giving ru-tum — 3 orbits for our fiducial parameters. We 
therefore have the ordering ru-tum < Tdiff < Tnb/2, which 
ensures that the hydrodynamical corotation torque remains 
close to its maximum fully unsaturated value. 



2.4 Numerical method 

We use tw o different MHD codes : the finite v olume code 
RAMSES l|Tevssierj 120021 : iFromang et~ail 120061 ) using the 
MUSCL-Hancok Godunov Scheme, and the finite difference 
code NIRVANA (jZiegler fc Yorkd ll997Tl which uses an a l- 
gorithm similar to the ZEUS code I Stone fc Normanll 19921 ). 
The results of the two codes are shown to be in good agree- 
ment in Appendix A. Figures [TJ to \5\ [7] [8] and [T2l to [TBI were 
done using the code RAMSES. Figures [6] and [9] were done 
using NIRVANA. Finally, Figures [9] [19] and [20] were done 
using both codes. 

The numerical domain extends around the planet loca- 
tion in the range r G [0.5,2], and (f> G [— 7r,7r]. Wave-killing 
zones are used near the grid's inner edge (r £ [0.5, 0.65]) and 
outer edge (r G [1.7, 2]) in order to avoid spurious reflections 
of the planet's wakes. For this purpose the density and ve- 
locity are damped to their initial values. The magnetic field 
however is not damped, because the constraint that its di- 
vergence remains zero would require a special treatment. We 
found that this was enough to avoid any significant refection 
at the edges where reflective boundary conditions are ap- 
plied (see Figure [T}. Our default grid resolution is n r = 512, 
714, = 1024 with RAMSES, and n r = 468, = 1248 with 
NIRVANA. It is such that the half-width of the horseshoe 
region is resolved by about 8 cells for our fiducial set of pa- 
rameters. In Appendix A, we show that this resolution is 
good enough to obtain converged results. 



3 MAGNETIC FIELD DISTRIBUTION INSIDE 
THE PLANET'S HORSESHOE REGION 

A global view of a simulation of Model 1 with the fiducial 
parameters is represented in Figure [T] The density pertur- 
bation (left panel) shows two distinct features: the usual 
wake (i.e. the spiral density waves launched in the inner and 
the outer disc), and an underdense region localised around 
the planet's orbital radius. This second feature is due to the 
horseshoe dynamics in the presence of a magnetic field and 
will be described further in Section f4.1l The magnetic field 
perturbation shows also two features related to the wake and 
the horseshoe dynamics, but the latter is clearly the most 
prominent one. This highlights the fact that the main effect 



of the relatively weak magnetic field considered here is on 
the horseshoe dynamics. 

We continue the description of the numerical results by 
considering the effect of the horseshoe gas motion on the 
magnetic field configuration, as is illustrated by Figures [2] 



for Model 1 (E 



-1/2 



, T oc r 1 ) , and by Figure [3] for 



Model 2 (£ oc r~ 3//2 , uniform temperature). We compare 
simulations using our fiducial magnetisation j3 = 100 (right 
panels) with simulations using an extremely weak magneti- 
sation (/? = 10 s , middle panels) where the magnetic field 
is effectively passive. This allows to disentangle the mag- 
netic field perturbation due to a pure advection/diffusion of 
the magnetic field from its dynamical effects on the velocity 
field. The extent of the horseshoe region is shown by the 
separatrices overplotted with white lines: as expected the 
magnetic field is not strong enough to prevent horseshoe 
motion of the gas. The magnetic configuration (illustrated 
by the magnetic pressure with the colour contours and the 
magnetic field lines in black) is quite similar in the passive 
and non-passive field cases: the magnetic field lines are bent 
by the horseshoe motion of the gas, resulting in an accumu- 
lation of magnetic flux along the downstream separatrices 
and in particular close to the stagnation point. The width 
of the region over which the magnetic flux accumulates is set 
by an equilibrium between the advection due to horseshoe 
motion bringing magnetic flux toward the downstream sep- 
aratrix and the diffusion of this magnetic flux by the action 
of resistivity. Decreasing the resistivity results in a narrower 
region of accumulation, and therefore a stronger amplifica- 
tion of the magnetic field there (not shown in the figures 
because the pattern is very similar). We note that the pres- 
ence of a non- negligible magnetic field induces slightly more 
compression of the magnetic field near the planet. This may 
be linked to a (rather surprising) slight increase of the width 
of the horseshoe region due to the presence of the magnetic 
field (see Figure I16[) . Indeed a wider horseshoe region re- 
sults in a larger shear velocity at the horseshoe separatrix 
and a larger advection velocity in the horseshoe region, mak- 
ing the advection process (and therefore the magnetic field 
compression) more efficient. 

Another point which will be of great importance in our 
discussion of the torque in Section [4] is the asymmetry of 
the configuration with respect to the planet azimuth. Fig- 
ure [2] and [3] show that the stagnation point is shifted in 
the azimuthal direction with respect to the planet location: 
in Model 1 the stagnation point has a negative azimuth, 
while it has a positive azimuth in Model 2. The sign of this 
asymmetry is the same in the MHD and the hydrodynam- 
ical cases, but the amplitude of the asymmetry is larger 
in the presence of the magnetic field. In Section [4] we will 
show by varying the disc profiles that this is a general re- 
sult. The origin of the stagnation point asymmetry in the 
hydrodynamic case will be further discussed in the next few 
paragraphs. For now, we note that the asymmetry in the 
stagnation point azimuth induces an asymmetry in the dis- 
tribution of magnetic pressure. Indeed the magnetic field is 
more concentrated near the stagnation point, and the mag- 
netic pressure is therefore maximum on the same side of the 
planet. 

Finally, we compare the MHD simulations with a pas- 
sive magnetic field (middle panels of Figures [2] and [3} to 
a passive evolution found by solving an advection/diffusion 
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Figure 1. Global view of a simulation of Model 1 (p = —1/2, q = —1) with the fiducial parameters. The left panel shows the surface 
density perturbation, while the right panel shows the perturbation of the azimuthal component of the magnetic field. The location of 
the planet is shown with a white + sign at r = 1, ip = 0. The extent of the region around the planet represented in Figures [21 E] and [5] is 
shown with a white sector in the left panel. 



kinematic calculation 





Figure 2. Magnetic field configuration in model 1 (p = —1/2, q = —1). Left panel: linear hydrodynamical calculation coupled with an 
advection-diffusion equation for a passive magnetic field. Middle panel: MHD simulation with a very weak magnetic field (/? = 10 s ). Right 
panel: MHD simulation with our fiducial magnetic field strength (/3 = 100). The colour contours show the magnetic energy normalised 
by its initial value. The black lines represent magnetic field lines, while the white lines represent the separatrices delimiting the horseshoe 
region. The black cross shows the position of the stagnation point at the intersection of the two separatrices. The planet location is 
depicted by a black + sign at r = l,tp = 0. 



equation for the magnetic field with a fixed velocity field 
(left panels) obtained from a linear calculation of the disc 
response to forcing by the planet. To perform these calcu- 
lations the velocity field was constructed from linear re- 
sponse calculations for azimuthal mode n umbers m up to 
60 (see iPaardekooper fc Papaloizou! 120081 ). Singularities in 
the linear response were dealt with using the Landau pre- 
scription which adds a purely imaginary frequency shift to 



the forcing frequency. The magnitude of the shift was taken 
to be 10~ 3 of the frequency itself, though as only a loga- 
rithmic singularity is involved in the worst case, changing 
this to 10 -5 was found to make no difference. This fea- 
ture affects the solution only on the corotation circle where 
the streamlines turn. This part of the flow requires a non- 
linear a nalysis, as it is not correctly repr esented by linear 
theory (|Paardekooper fc Papaloizoull200a) . Nonetheless, al- 
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Figure 3. Same as Figure[2]but for model 2 (p = —3/2, q = 0). 



though the horseshoe region is somewhat narrower than that 
found from nonlinear simulations, as can be seen from Fig- 
ures [5] and the procedure reproduces the correct asymme- 
try in the location of the stagnation point. Such an asym- 
m etry was already observe d in hydrodynamical simulations 
by Casoli & Massct (2009) and Paardckoopcr & Papaloizou 
(200$. iPaardekooper fc Papaloizoul (|2009h showed that it 
was due to an asymmetry between the inner and outer 
wakes. Indeed they showed that the azimuthal asymmetry of 
the stagnation point disappears if the gravity of the planet 
is cut off at distances larger than H, which prevents the ex- 
citation of the wake while keeping the horseshoe dynamics. 
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In order to determine the kinematic effect on an initially 
purely azimuthal field, we solved the induction equation us- 
ing the velocity field calculated with the linear response of 
the disc to the planet, and with the same magnetic diffusiv- 
ity as in the full MHD simulations. The numerical domain 
is centred on the planet and extends over the full 2ir in the 
azimuthal direction, and over a length 0.5r p in the radial 
direction. Other parameters of the calculations were as in 
models 1 and 2. For the problem on hand the induction 
equation reduces to an advection diffusion equation for the 
flux function that is readily soluble as an initial value prob- 
lem. At azimuthal boundaries we imposed periodic bound- 
ary conditions and at radial boundaries the magnetic flux 
was held fixed at its initial value, so conserving the total 
azimuthal magnetic flux. The ultimately steady state solu- 
tions illustrated in Figures [2] and [3] show good qualitative 
and quantitative agreement with the non linear simulations 
for weak fields in that the field gets concentrated along out- 
going separatrices with an azimuthal asymmetry in the field 
strength that produces a larger field on the side of the planet 
containing the stagnation point. The important feature is 
the asymmetry of the stagnation point in azimuth and the 
associated asymmetry in field strength, which seems to be 
the main determinant of the corotation torque excess (see 
Section 3]). This is quite well reproduced by linear calcula- 
tions and being a linear feature, it should apply generally to 
the low mass planet regime. 
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10 20 bO 40 bO 60 

time In orbits 

Figure 4. Time evolution of the torque exerted on the planet in 
models 1 and 2 (shown in black and red respectively), with and 
without a magnetic field (full and dashed lines respectively). 



4 AN ADDITIONAL COROTATION TORQUE 
OF MAGNETIC ORIGIN 

4.1 Two particular cases 

We now describe the torque exerted by the disc on the planet 
in the MHD and corresponding hydrodynamical simulations, 
and start by considering the two fiducial models described in 
Section [3] As shown in Figure |4j the torque is dramatically 
changed in the presence of a magnetic field. An additional 
torque component is observed in both Model 1 and 2, which 
takes about 10-15 orbits to establish. This corresponds to a 
few times the horseshoe U-turn timescale, which is about 3 
orbits. This additional corotation torque has a large value 
with these fiducial parameters, and the total torque may 
become positive (in model 1). This underlines the potential 
importance of a weak magnetic field on type I migration, 
since here the magnetic pressure is only one percent of the 
thermal pressure. 

While the amplitude of the additional torque is very 
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Figure 5. Comparison of the density perturbation (i.e. the density in the stationary state minus the density in the initial state) in 
hydrodynamical (left panels) and MHD simulations (right panels). Results are shown for model 1 (upper row), and model 2 (lower row). 
The black lines represent the streamlines, the white lines the separatrices of the horseshoe region, and the black cross the stagnation 
point. 
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similar in both models, its sign differs: it is positive in 
Model 1 and negative in Model 2. To investigate this sign 
difference, we show the density perturbation induced by the 
planet for both models in Figure [5] Comparing the hydro- 
dynamical simulations (left panels) and MHD counterparts 
(right panels) , we observe that the magnetic field induces un- 
derdense lobes located around the downstream separatrices. 
This can be interpreted as being due to the accumulation 
of magnetic pressure at the same location, which leads to 
a decrease of thermal pressure (and thus of surface density) 
to maintain approximate total pressure balance. We also no- 
tice a clear rear/front asymmetry in the density perturba- 
tion inside the planet's horseshoe region. The amplitude of 
the negative density perturbation takes its maximum value 
on the side of the planet where the stagnation point lies: 
at negative azimuth for Model 1, and positive azimuth for 
Model 2. This is consistent with the perturbation of mag- 
netic pressure being maximum about the stagnation point. 
The origin and the sign of the additional torque can then 
be interpreted with this asymmetry of the density distribu- 
tion. An underdense region induces a negative torque if it 
is in front of the planet (i.e. at positive azimuth), and a 
positive torque if it lies behind the planet (i.e. at negative 
azimuth). When the stagnation point has a positive azimuth 
as in Model 2, the largest underdensity lies in front of the 
planet and the additional torque is therefore negative. Con- 
versely, when the stagnation point has a negative azimuth 
as in Model 1, the largest underdensity is behind the planet 
and the additional torque is therefore positive. 

Although the above interpretation of the MHD torque 
excess explains well the sign of the torque, we should add 
a note of cau tion in interpreti n g the torque using density 
lobes. Indeed, iMasset fc Casola (|2009l ) showed that such an 
interpretation could be misleading in the case of the corota- 
tion torque in an adiabatic hydrodynamical disc, because a 
diffuse density perturbation plays an important role as well. 
Instead they argued that the torque should be interpreted 
as arising from a creation of vortensity at the separatrices of 
the horseshoe region. In the case we are considering here, the 
magnetic field is also responsible for a creation of vortensity. 
This might be used to give an alternative explanation of the 
torque excess, although it is not yet clear whether a similar 
interpretation applies in the presence of a magnetic field. 
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4.2 Dependence on the gradients of density, 
temperature and magnetic field strength 

We have shown in the previous section that a new compo- 
nent of the corotation torque arises in the presence of a weak 
magnetic field. Its amplitude is the same but its sign differs 
in two disc models that differ by the choice of the temper- 
ature and density profiles. To better characterise this new 
torque component, we performed series of MHD and HD 
runs in which we varied the slopes of the background mag- 
netic field, density and temperature profiles. The results of 
these series of run are presented in the next subsections. 



4-2.1 Varying the magnetic field strength gradient 

Strikingly, the torque does not depend at all on the gradient 
of magnetic field strength ! This is illustrated in Figure [(J 



Figure 6. Results of simulations with different slopes of the ini- 
tial magnetic field profile : b = — 1 (black), b = —0.5 (blue) and 
6 = (red). Upper panel : azimuthally averaged magnetic field 
profile in a steady state normalised by its initial value at the 
planet location (the vertical dashed lines show the extent of the 
horseshoe region). Lower panel : time evolution of the torque. 
Strikingly the torques from the different simulations are indistin- 
guishable. 



showing the azimuthally averaged magnetic field strength 
(upper panel) and the time evolution of the total torque 
(lower panel). Three models with different slope of the mag- 
netic field strength are shown : b = —1,-0.5 and 0. The 
magnetic field profiles differ substantially (except within 
the planet's horseshoe region) but the time evolution of the 
torque in the different models is indistinguishable. 

This surprising result contrasts with the usual be- 
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Figure 7. Time variation of the MHD torque excess (defined as 
the torque in the MHD simulation minus the torque in its hy- 
drodynamical counterpart) for various values of the temperature 
slope index <j, and a density profile p = —1/2. 



haviour of the corotation torque, which is often proportional 
to the gradient of some conserved quantity (vortensity in an 
isothermal disc, entropy for a disc with an energy equation). 
It also differs from the regime of strong mag netic field where 
the torque arises from magnetic resonances : iTeraueml (|2003h 
indeed found that the sign and magnitude of the torque de- 
pends strongly on the slope of the magnetic field strength 
profile. How should this different behaviour be interpreted 
? One might interpret the independence on the magnetic 
field gradient by the fact that the horseshoe motion redis- 
tribute the magnetic flux in the horseshoe region (it tends 
to accumulate it on the downstream separatrix). Therefore 
what matters is the total flux in the horseshoe region rather 
how it is initially distributed in it. The usual proportion- 
ality to gradients is due to the fact that gradients are the 
cause of a symmetry breaking, necessary for the presence of 
a net torque. Here the presence of a magnetic field induces 
an azimuthal asymmetry of the horseshoe motion and of the 
density distribution, which is not proportional to any gra- 
dient. Indeed, it will be shown later that the sign of this 
asymmetry is governed by the gradients of density and tem- 
perature, but its amplitude depends only on the magnetic 
field strength. Thus the MHD torque excess is not propor- 
tional to any gradient because an asymmetry seems to be 
caused by the interaction between the magnetic field and 
the horseshoe motion. 
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4-2.2 Varying the temperature gradient 

Figure [7] shows the time evolution of the MHD torque ex- 
cess (defined as the torque in the MHD simulation minus 
the torque in its hydrodynamical counterpart) for a series 
of runs varying the temperature gradient. It shows that the 
torque excess changes sign abruptly around q = 0, while its 
amplitude in steady state remains roughly independent of 
q within 10-20%. This behaviour is also displayed on the 
top and middle panels of Figure [H] showing respectively the 
torque and the MHD torque excess as a function of the index 
of the temperature profile q. Importantly, this discontinuous 



Figure 8. Results of a series of simulations varying q for p = 
— 1/2. Upper panel: total torque in a steady-state as a function 
of the slope of the background temperature profile, for MHD (red 
stars) and HD simulations (black diamonds). Middle panel: MHD 
torque excess, defined as the MHD torque subtracted by the hy- 
drodynamical torque. Red star symbols show the absolute value 
of the torque excess (when it is negative). Lower panel: Azimuth 
of the stagnation point in hydrodynamical (black) and MHD sim- 
ulations (red). The dependence of <p B with the temperature gra- 
dient is simil ar to that of the outer X stagnation point obtained 
in Fig. 10 of ICasoli fc Massed j200ijh (the other two stagnation 
points appear only at smaller softening length, as discussed in 
Section I4.4|l . 
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Figure 9. Same as the middle panel of Figure \E\ but with p = 
—3/2: MHD torque excess (defined as the MHD torque subtracted 
by the hydrodynamical torque) as a function of the temperature 
slope index q. Red star symbols show the absolute value of the 
MHD torque excess (when it is negative). For comparison the 
dotted line shows the hydrodynamical torque. 



behaviour of the torque excess follows that of the azimuth 
of the stagnation point, shown in the lower panel of Fig- 
ure [8] Indeed the azimuth of the stagnation point changes 
sign abruptly around the same temperature gradient as the 
torque excess (q = 0), in the MHD as well as the hydro- 
dynamical simulations. This confirms the interpretation of 
the MHD torque excess given in Section 14.11 as an asym- 
metry of the density distribution in the horseshoe region 
driven by the azimuthal shift of the stagnation point with 
respect to the planet. It also confirms that the sign of the 
stagnation point azimuth is the same in the MHD and the 
hydrodynamical cases, as was noted in the particular cases 
of Model 1 and 2. The amplitude of the azimuthal shift is 
however different: in the MHD case it appears to be indepen- 
dent of the temperature gradient, while in the hydrodynami- 
cal case it decreases when approaching the transition where 
it changes sign. Thus the sign of the stagnation point az- 
imuth is governed by a hydrodynamical process determined 
by the density and temperature gradients (linked with an 
asymmetry between the inner and outer wakes as discussed 
in Section [3|, while its amplitude seems to be driven by a 
magnetic process independent of these gradients (its depen- 
dence on the magnetic field strength will be discussed in 
Section T4.3p . Finally, we note that while the final amplitude 
of the torque barely depends on the temperature gradient, 
the time evolution to reach this amplitude differs signifi- 
cantly: the closer to the transition where the torque excess 
changes sign abruptly, the longer it takes for the torque ex- 
cess to build up (Figure [7|). This may be interpreted in the 
following manner. The presence of the magnetic field tends 
to amplify any azimuthal asymmetry present in the hydro- 
dynamical flow to finally reach an amplitude independent of 
the initial asymmetry. It may therefore seem natural that a 
smaller initial asymmetry takes longer to be amplified by the 
magnetic field, resulting in the time evolution of the torque 
excess observed in Figure [7] 

Figure [9] shows the dependence of the MHD torque ex- 
cess with the temperature gradient using a different density 
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Figure 10. Results of a series of runs where the density profile 
is varied at fixed uniform temperature profile (q = 0). The MHD 
torque excess is shown as a function of the slope index of the den- 
sity profile p. Red symbols show the absolute value of the torque 
excess. For comparison the dotted line shows the hydrodynami- 
cal torque. The simulations were performed using both the codes 
RAMSES (diamonds) and NIRVANA (plus signs), which allows 
to assess the numerical errors. The two codes agree within 5% for 
the amplitude of the torque excess, and find the sign transition 
at q= -0.6 ± 0.2. 



profile (p — —3/2). This confirms the qualitative features 
observed with p = —1/2: the amplitude of the torque ex- 
cess is roughly constant, but its sign changes abruptly at a 
critical temperature gradient. The value of the temperature 
gradient at the transition depends on the density gradient 
(it is q ~ -0.55 for p = -3/2, and q ~ -0.05 for p = -1/2) 
but it still corresponds to the transition where the stag- 
nation point azimuth changes sign in both hydrodynamical 
and MHD simulations (not shown here). The dependence 
on the density profile is further investigated in the following 
subsection. 



4-2.3 Varying the density gradient 

We performed a series of runs varying the density profile at 
fixed uniform temperature profile (q — 0). Figure [TU] shows 
the dependence of the MHD torque excess with the slope 
index of the density profile p. The behaviour is very similar 
to that described in Section [4.2.21 when varying the temper- 
ature profile: the amplitude of the torque excess is barely 
dependent on the density profile, but its sign does depend 
on it. A sharp transition from a negative to a positive MHD 
torque excess is obtained at p — —0.55 ±0.15 (between —0.7 
and -0.6 for NIRVANA, and between 0.5 and 0.4 for RAM- 
SES). Therefore, steeply decreasing density profiles give a 
negative MHD torque excess, while shallower (or increasing) 
profiles drive a positive MHD torque excess. 

Figure [TU] can also be used to assess the importance of 
numerical errors, since the series of simulations has been per- 
formed both with the RAMSES and NIRVANA codes. We 
find that the qualitative behaviour of constant amplitude 
and sharp sign transition is a robust feature of the simula- 
tions. The amplitude of the torque excess agrees within 5% 
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Figure 12. Diagram showing stable (blue diamonds) and un- 
stable (red diamonds) runs in the {a, /3} plane. The dashed line 
corresponds to a fit of the marginal stability: ct\fp = 3.10 - 2 . 



Figure 11. Azimuth of the (outer-X) stagnation point (<p s ) in a 
series of non-magnetised simulations, as a function of the back- 
ground surface density and temperature gradients at the planet 
orbit. The solid line shows the linear function p — 1.3q + 0.6 = 0. 
Wherever tp 3 < 0, the magnetic-related horseshoe drag is positive, 
while it is negative for ip s > 0. The position of the two fiducial 
models is represented with a white circle. The white rectangle 
represents the range of parameters suggested by observations of 
protoplanetary discs at large radii: q 6 [—0.7,-0.4], p G [—1,0] 
(see Section 16.21 for a discussion). 



between the two codes, and the location of the sign transi- 
tion agrees within Ap = 0.2. 



4-2.4 Sign of the MHD torque excess 

We have found that the amplitude of the MHD torque ex- 
cess is roughly independent of the temperature, density and 
magnetic field strength gradients. On the other hand, the 
sign of the torque excess does depend on both the tempera- 
ture and density gradients. This sign is dictated by the sign 
of the stagnation point azimuth (which is the same in both 
hydrodynamical and MHD simulations): a positive stagna- 
tion point azimuth leads to a negative torque excess, while 
a negative azimuth leads to a positive torque excess. This is 
interpreted as due to a rear-front asymmetry of the under- 
dense lobes near the downstream separatrices, the strongest 
underdensity arising on the side of the planet where the 
stagnation point lies. 

Given that the azimuth of the stagnation point is ob- 
served to have the same sign in HD and MHD simulations, 
one can use HD simulations to infer its sign and thus the sign 
of the MHD torque excess. For this purpose, we have run 
HD simulations varying both the density and temperature 
profiles. The simulations were performed with FARGO: a 
two-dimensional hydrodynamical code similar to NIRVANA, 
which uses the so-called Fargo algorith m to speed u p the 
calculation by increasing the time-step ljMassetH200(jh . The 
results are presented in Figure [TT] showing the azimuth of 



the stagnation in the plane {p, q}. A region of positive az- 
imuth at low p and large q is clearly separated from a region 
where the azimuth is negative at large p and low q. The limit 
between these two regions of parameter space can be approx- 
imated by a linear function of p and q: p — 1.3g + 0.6 = 0. 
Thus we find that the sign of the stagnation point azimuth 
is given by: 



sign(<£ a ) = sign(1.3g - p - 0.6). 



(17) 



As mentioned before, the MHD torque excess has the oppo- 
site sign. This approximate formula can now be compared 
to the results of Section |4T] |4~2~21 and 14.2.31 We first note 
that Model 1 (respectively 2) lies in the region of negative 
(positive) stagnation point azimuth, in agreement with the 
results of Section 14.11 At p = —0.5, the formula predicts a 
sign transition at q = 0.08, to be compared with q ~ —0.05 
found in Section |4~2~21 ( Figure l8jl . At p = —3/2, the formula 
gives q = —0.7 against q ~ —0.55 for the simulations of Fig- 
ure[5] Finally, at q — the formula predicts a sign transition 
at p — —0.6 against p = —0.55 ± 0.15 for the simulations of 
Figure 1101 All in all, equation \T7\ therefore reproduces the 
results of the simulations with a precision of 0.1 — 0.2 in p 
and q. 



Dependence on magnetic field strength and 
diffusion coefficients 



4.3 



As shown in Section [4.21 the amplitude of the MHD torque 
excess does not depend on the gradients of magnetic field 
strength, density and temperature. It may therefore depend 
only on the magnetic field strength (through /3), the diffusion 
coefficients (a and V), the mass of the planet (M p /M»), the 
aspect ratio of the disc (h) and the softening length of the 
planet potential. In this section we explore the dependence 
on the magnetic field strength and the diffusion coefficients, 
and leave for future work the study of different values of the 
planet mass and aspect ratio of the disc. The effect of varying 
the softening length is briefly described in Section \4. 41 
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Figure 13. Amplitude of the MHD torque excess as a function of 
P for three different values of a: 2.10 - 3 (red), 5.10 - 3 (black) and 
10 -2 (green). The dashed lines show the fit given by Equation ll9l 
(which assumes V oc l/(/3 2 cv 4 )). For strong magnetic fields or low 
diffusion coefficients the torque excess deviates from this scaling. 
Other parameters are the fiducial ones and the disc Model 2. 
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Figure 15. Amplitude of the MHD torque excess as a function of 
the magnetic Prandtl number. The black diamonds show a series 
of run where the viscosity is kept constant to a = 5.10 - 3 , while 
the resistivity is varied. Reversely, the red squares correspond to 
the resistivity being constant a/V = 5.10 -3 , while the viscosity 
is varied. The magnetic torque excess has a steep dependence on 
the resistivity, but only a very weak one on the viscosity (for the 
range of diffusion coefficients considered). Other parameters are 
the fiducial ones and the disc Model 2. 
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Figure 14. Amplitude of the MHD torque excess as a function 
of a for p = 100. Other parameters are the fiducial ones and the 
disc Model 2 (black diamonds) or Model 1 (red + signs). Here, 
the dimensionless resistivity a/V (which we show in Figurc [T5l to 
be the relevant diffusion parameter determining the MHD torque 
excess) is equal to a with our fiducial value of V = 1. The dashed 
line shows the fit given by Equation ll9l with T oc (V/a) 4 . For low 
values of a the torque excess deviates from this scaling. 



4-3.1 Stability domain 

When varying a and /3, we observed that for small values of 
these parameters the horseshoe motion could become unsta- 
ble. This instability grows on the downstream separatrices 
of the horseshoe region. It is therefore due to the presence of 
the planet rather than an intrinsic property of the disc. As 
an additional check, we ran simulations without a planet but 
with small velocity perturbations added to the background 
flow, and found no instability. The instability is probably 
caused by a creation of vortensity due to the interaction of 
the magnetic field with the horseshoe motion. When the 



vortensity perturbation is strong enough, this leads to a 
shear instability that creates vortices near the downstream 
separatrices. A detailed description of this instability is how- 
ever postponed to future work. Here we focus on the stable 
region of parameter space to characterise the dependence of 
the torque excess on both a and /3. 

The region of parameter space {a, j3} where the horse- 
shoe motions are stable or unstable is shown in Figure 1121 
From this figure, we deduce the following criterion for the 
motion to be stable for our fiducial values of planet mass 
and aspect ratio: 



xsj~fi > 3 x 10" 



(18) 



In the following Sections, we present the results of simulation 
in this parameter domain. 



4-3.2 MHD torque excess 

Figure [13] shows the amplitude of the MHD torque excess 
in three series of runs varying /3 for three different values 
of a. As expected for a magnetic effect, the MHD torque 
excess increases with the magnetic field strength (i.e. when 
/3 decreases). It also increases very steeply when the diffu- 
sion is decreased, as can be seen by comparing the different 
series in Figure [13] and more directly in Figure [14] showing 
the dependence of the MHD torque excess on a for ft = 100. 
The latter also shows that the MHD torque excess vanishes 
for high values of the diffusion coefficients, as can be ex- 
pected because the magnetic field and gas motion decou- 
ple when the resistivity is very large. The dependence on 
the diffusion can be interpreted in the following way. The 
MHD torque excess is due to an azimuthally shifted accu- 
mulation of magnetic pressure along the downstream sepa- 
ratrix. The efficiency with which the magnetic flux can be 
accumulated at the separatrix is governed by an equilibrium 
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between the advection by the horseshoe motion and the dif- 
fusion by the resistivity: the lower the resistivity the thinner 
the layer where the magnetic flux is accumulated and there- 
fore the stronger the magnetic pressure there. This stronger 
magnetic pressure leads to a stronger underdensity, thus ex- 
plaining the steep increase of the MHD torque excess when 
the diffusion is decreased. 

Following this reasoning, we expect that the steep de- 
pendence of the torque excess is mostly due to the resistivity 
rather than the viscosity (both of which are varied concur- 
rently in these runs where the magnetic Prandtl number is 
fixed to V = 1). We checked this by performing two se- 
ries of runs where one of the diffusion coefficients is varied 
while keeping the other one constant and show the results 
in Figure [TJ] as a function of the magnetic Prandtl num- 
ber. Varying the viscosity by a factor 10 changes the MHD 
torque excess by no more than 20%, if the resistivity is kept 
constant. By contrast, a similar variation of the resistivity 
at fixed viscosity changes the torque excess by a factor 40 ! 
To first order, the MHD torque can therefore be considered 
to be independent of the viscosity in the range of parameters 
explored in this study. This is probably due to the fact that 
the values of the viscosity considered are close to that lead- 
ing to a fully unsaturated horseshoe drag, around which the 
dependence on the viscosity is quite weak. It is expected that 
decreasing the viscosity further would eventually lead to a 
saturation of the corotation torque as in the hydrodynami- 
cal situation, due to a lack of angular momentum exchange 
between the coorbital region and the rest of the disc. 

Figure [13] suggests that the MHD torque excess has a 
power law dependence with respect to /3 for moderately large 
values of j3 and a. We find that the torque excess is fairly 
well fitted with the following scaling (for our fiducial values 
of the planet mass and disc aspect ratio): 
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Figure 16. Characteristics of horseshoe region as a function of 
for a = 5.10~ 3 and model 2. Upper panel: azimuth of stagnation 
point. The dotted line shows the value for a non magnetised run. 
Lower panel: half width of horseshoe region measured at an az- 
imuth of ip = ±1. The rear (ip = —1) is shown in black, while the 
front (ip = 1) is shown in red. Other parameters are the fiducial 
ones and the disc Model 2. 



This fit is overplotted in Figures [13] and [14] with dashed 
lines. It is valid for large enough values of a and /3 corre- 
sponding to a MHD torque excess smaller or comparable to 
the torque in a hydrodynamical case. For smaller values of 
these two parameters, the dependence flattens out leading 
to a weaker torque excess than predicted by the above scal- 
ing. It is interesting to note that this scaling depends on the 
same combination of powers of a and /3 as in the stability 
criterion (equation 118}) : the value of the dimensionless pa- 
rameter /3a 2 determines both the stability of the horseshoe 
motion and the amplitude of the MHD torque excess. Al- 
though the physical origin of this scaling remains somewhat 
uncertain and should be investigated further in the future, 
we propose the following tentative interpretation. The den- 
sity perturbation in the horseshoe region can be expected 
to be inversely proportional to /3 as it is proportional to 
the magnetic pressure. The torque excess scaling as 1//3 2 
may then be understood if the rear-front asymmetry of the 
density lobes is also inversely proportional to /3 (see below 
the dependence of the stagnation point azimuth with /3). 
One may further speculate that the density perturbation 
and rear- front asymmetry scale as 1/a 2 , leading to the de- 
sired scaling of the torque excess. 



4-3.3 Horseshoe region 

After studying the torque, we turn to the characteristics of 
the horseshoe region and their dependence on the magnetic 
field strength. Fieure [T6l shows the azimuth of the stagnation 
point and the half-width of the horseshoe region as a func- 
tion of /3. At large /3, the azimuth of the stagnation point and 
the half width of the horseshoe region tend to their hydrody- 
namical value, as expected in the weak magnetic field limit. 
As the magnetic field increases, the azimuth of the stagna- 
tion point increases significantly, confirming the idea that 
the presence of the magnetic field amplifies the azimuthal 
asymmetry of the horseshoe motion. The width of the horse- 
shoe region also increases with stronger magnetic fields (in 
this weak field regime). This result is somewhat surprising 
because one may naively expect the magnetic tension to re- 
sist the horseshoe motion, and therefore that the width of 
the horseshoe region decreases with stronger magnetic field. 
The reason for the observed behaviour remains unclear. We 
also observe that the rear-front asymmetry of the horseshoe 
width is more and more pronounced as the magnetic field is 
increased. 
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Figure 17. MHD torque excess (diamonds) and hydrodynam- 
ical torque (stars) as a function of the softening length of the 
planet's gravitational potential. The simulations were performed 
using Model 2 and the fiducial parameters. Red symbols cor- 
respond to torques computed excluding the Hill sphere of the 
planet, while it is not excluded in the torques shown by the black 
symbols. The dashed and dash-dotted lines show a dependence 
as 1/e for comparison. The vertical dotted line shows the value of 
the planet's Hill radius (given by equation \12) . 



4.4 Dependence on the softening length 

So far all the results presented in this paper have used a 
softening e — 0.6H(r p ) for the gravitational potential of the 
planet. It is not clear which value of the softening length is 
the best representation of the 3D behaviour in discs, there- 
fore it is useful to study the dependence of our results with 
this parameter. Figure [17] shows the dependence of the hy- 
drodynamical torque and of the MHD torque excess on the 
softening length. As the softening length is decreased, both 
torques increase significantly. The MHD torque excess has a 
dependence that is roughly consistent with being inversely 
proportional to e in this range of softening length, compara- 
ble to the scaling of the hydrodynamical corotation torque in 
locally isothermal discs (jPaardekooper et al.ll201ol ). We also 
considered the effect of excluding the planet's Hill sphere 
from the torque calculation (which was not done so far in 
this paper). This has a significant impact only if the soft- 
ening length is smaller than the Hill radius, and therefore 
barely affects the results for our fiducial softening length. 

Another point of interest is the behaviour of the stag- 
nation point as the softening length is decreased. Indeed, 
for small enough values of the softening length there may 
actua lly be three stagnation points in the hydrodynamical 
case (|Casoli fc M assct 20091). Given the importance of the 
stagnation point in our interpretation of the MHD torque 
excess, we may wonder how such a behaviour may affect our 
results. In the case of a softening length of e — 0.3H, we do 
observe the appearance of three stagnation points in our hy- 
drodynamical simulatio ns in agreement with the results of 
ICasoli fc Massed j2009h . However, the corresponding MHD 
simulation still shows only one stagnation point and a very 
similar topology as the simulations with a larger softening 
length. The sign of the stagnation p oint azimuth is th e n tha t 
of the outer X-point (as defined bv lCasoli fc Massed ()2009l ) 
as the intersection of the outer separatrices, which delimit 



the horseshoe region). The MHD dynamics therefore seems 
to depend on the azimuthal asymmetry of the hydrodynam- 
ical situation but not on the details of the topology of the 
streamlines. 



5 COMPARISON WITH 3D MHD 
SIMULATIONS OF DISCS WITH 
TURBULENCE DUE TO THE MRI 

In this section we compare our results with the 3D MHD 
simulations of discs performed by BFNM11, where the MRI 
operates throughout the whole disc. The goal is to check 
how well a 2D description where the effects of turbulence 
are modelled by diffusion coefficients can explain a 3D tur- 
bulent situation. For comparison purposes, we ran 2D simu- 
lations using the same parameters as in BFNM11. We study 
two disc models differing in the slope of the temperature and 
surface density profiles: model 1 (p — —1/2 and q = — 1) and 
model 2 (p = —3/2 and q — 0) as defined earlier in this pa- 
per. The other parameters are: M p — 3x 10 -4 , H/r = 0.1, as 
well as ft = 50 and a = 3x 10~ 2 as measured in the turbulent 
simulations. The simulations of BFNM11 are not vertically 
stratified and use a softening length of the planet potential 
e = 0.2H(r p ), which we adopt here. With this rather small 
value of the softening length, the Hill sphere can have an 
effect on the torque (see Section [4.4p . For consistency with 
BFNM11, we exclude the Hill sphere from the torque cal- 
culations. In a 3D stratified simulation, the structure of the 
flow may differ somewhat from that obtained with a soft- 
ened potential. This should be considered in the future for a 
better quantitative estimate of the MHD torque excess. The 
slope of the magnetic field strength, which was found to have 
no effect on the MHD torque excess (see Section [4727T]) , was 
given our fiducial value 6 = — 1, which avoids diffusion since 
it is current-free. The effective magnetic Prandtl number 
was not measured in BFNM11 because it is not a straight- 
forward analysis and requires to run specific simulations. As 
a result, the magnetic Prandtl number needs to be specified 
based on expectations from the literature. As discussed in 
Section \2. 11 it is expected to be of order unity but may differ 
from it by a factor of a few. We therefore considered three 
values of this parameter in order to determine its influence 
on the results : V = 2, V = 1, and V = 0.5. In order to 
avoid a viscous evolution of the surface density profile, we 
used a radial profile of the viscosity such that the radial ve- 
locity vanishes (a uniform value of the diffusion coefficients 
in model 1, and an a parameter scaling like r- 1 ' 2 in model 
2). In addition to the MHD simulations, we also ran hydro- 
dynamical simulations with the same parameters. 

Figure[TS]illustrates the overall flow topology of Model 1 
in turbulent and laminar simulations, by showing the density 
perturbation and the streamlines. The width of the horse- 
shoe region is in good agreement between the turbulent and 
the laminar simulation. In both simulations, we clearly ob- 
serve an underdense lobe near the downstream separatrix 
at negative azimuth. A much weaker underdensity is also 
marginally seen at positive azimuth. We showed in Section[3] 
and 14. II that these underdense lobes are due to an accumu- 
lation of magnetic field near the downstream separatrix of 
the horseshoe region, and that the azimuthal asymmetry of 
this underdensity (clearly seen here) leads to a net torque 
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Figure 18. Surface density perturbation: comparison of the 3D turbulent simulations of BFNM11 (left panel) and our 2D laminar 
simulations (right panel). The 3D simulations have been averaged in the vertical direction and on a time-period of 35 orbits between 
520 and 555 orbits. The azimuthally averaged surface density has been subtracted to obtain the density perturbations, shown with the 
colour contours. The streamlines are shown with black lines, and the horseshoe separatrices are represented by red lines. We show the 
result of Model 1 (p = —1/2 and q = —1), and the 2D simulations use a magnetic Prandtl number of V = 0.5 (the result with V = 1 is 
qualitatively very similar with only a slightly larger underdensity near the horseshoe separatrix) . 



on the planet. Importantly, the stagnation point is located 
at a negative azimuth in both the turbulent and the laminar 
simulations, confirming that the most prominent underdense 
lobe lies on the same side of the planet as the stagnation 
point. This comparison shows that the overall flow topology 
and the density perturbation is in good agreement between 
the turbulent simulations of BFNM11 and our laminar sim- 
ulations. 

The results of the torque calculations are illustrated by 
Figure 1191 showing the time evolution of the torque (up- 
per row) and the torque density distribution (middle row 
for our simulations, lower row for the turbulent simulations 
of BFNM11). The range of values of the total torque in 
steady state obtained in the MHD turbulent simulations of 
BFNM11 is shown in the upper panels with a grey shaded 
area. Several features show a very nice qualitative agree- 
ment between our laminar simulations and the turbulent 
simulations of BFNM11. First the sign of the MHD torque 
excess: as earlier in this paper we find a positive torque ex- 
cess in Model 1 and a negative one in Model 2. This agrees 
with the results of BFNM11 who find a marginally positive 
torque excess in Model 1 and a negative one in Model 2. 
Secondly the shape of the torque density distribution shows 
two additional features in the presence of a magnetic field 
as compared to the hydrodynamical simulation: a negative 
contribution in the inner part of the horseshoe region and 
a positive one in the outer part, as was already observed 
by BFNM11. These are due respectively to the underdense 
lobes in front and behind the planet, which are located along 



the downstream separatrices where the magnetic flux accu- 
mulates (see Sections 131 and l47Tj) . We also observe that the 
two peaks are not symmetric about the planet radius or 
the corotation radius. They are shifted to smaller radii in 
Model 1 and to larger radii in Model 2, again in agreement 
with BFNM11. This asymmetry can be related to the az- 
imuthal shift of the stagnation point with respect to the 
planet, since the underdense region lies along the down- 
stream separatrices. In Model 1 the stagnation point az- 
imuth is negative, such that the underdense region lies be- 
hind the planet on a larger proportion of the horseshoe re- 
gion and therefore the positive torque contribution extends 
to radii smaller than the corotation radius. Conversely, in 
Model 2 the stagnation point azimuth is positive, the un- 
derdense lobe is preferentially in front of the planet, and the 
negative torque contribution extends to radii larger than 
the corotation radius. The qualitative agreement (in the 
torque excess, the torque density distribution, the surface 
density perturbations, and the azimuth of the stagnation 
point) gives confidence that the physical picture described in 
this paper applies to the turbulent simulations of BFNM11. 

Quantitatively, however, the agreement between the 
torque calculations in laminar and turbulent simulations is 
not so good, though the comparison is made difficult by the 
strong dependence on the magnetic Prandtl number. Start- 
ing the comparison with the fiducial value V = 1, we ob- 
serve that the amplitude of the MHD torque excess is too 
large in Model 1 and too small in Model 2 as compared to 
BFNM11. This discrepancy can be resolved by choosing the 
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Figure 19. Comparison of our 2D laminar simulations with the 3D turbulent simulations of Baruteau et al. (2011, BFNM11). Two 
disc models are considered: Model 1, i.e. p = —1/2 and q = —1 (left column, simulations performed with RAMSES), and Model 2, 
i.e. p = —3/2, q = (right column, simulations performed with NIRVANA). Upper row: total torque as a function of time. The range 
of steady state torque values obtained in the MHD turbulent simulations of BFNM11 is illustrated with a grey shaded area. Middle 
row: torque density distribution in the steady state in our laminar simulations. Lower row : torque density distribution obtained by 
BFNM11 with the code NIRVANA (the radius has been normalised by the orbital radius of the planet, and the torque density by the 
maximum of the hydrodynamical density distribution). In all panels hydrodynamic simulations are shown with dashed lines, while their 
MHD counterpart are shown with full lines. The viscosity and magnetic field strength are set to a = 3.10 -2 , and /3 = 50 as measured 
in BFNM11. Furthermore, we considered 2 values of the magnetic Prandtl number (setting the resistivity) : P m = 1 (black lines) and 
Pm = 0.5 (red lines). The torques and torque densities are calculated excluding the Hill radius for consistency with BFNM11. 
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right magnetic Prandtl number, however a different value is 
needed in both models: V < 0.5 in Model 1, and V — 1 — 2 
in Model 2. This is not entirely satisfactory at the present 
because we can give no reason for the magnetic Prandtl 
number of the MHD turbulence to be different in the two 
disc models. However, this may indicate a sensitivity of the 
results to the way turbulent diffusion operates and this may 
differ in the two cases. We note that the range of magnetic 
Prandtl numbers needed to explain the results of BFNM11 
is in line with numerical simulations of turbulence driven 
by the MRI, where the effective magnetic Prandtl number 
is measured to be in the range V ~ [0.3 — 5] dependin 



on the magnetic field configuration ( Guan fc Gammig |200 
iLesur fc Longarettill2009l ; iFromang fc Stonell2009l ). 



Comparing the laminar simulations presented in this 
section to the results of the previous sections, an unexpected 
feature appears: the amplitude of the MHD torque excess 
differs by a factor slightly more than 3 between the two disc 
models, while it was found in Section 14.21 that the MHD 
torque excess barely depends on the density and temper- 
ature profiles. This different behaviour is probably caused 
by non-linear effects, which are much more pronounced here 
due to the larger planet's mass (M p /(M»/i 3 ) = 0.3) and the 
smaller softening length. Indeed, we checked that increasing 
the softening length to e = 0.61/ reduced the discrepancy 
between the two models (a factor 2 in the MHD torque ex- 
cess). Also, we checked in Figure [T4l that the result of Sec- 
tion [472] was not a peculiarity of the fiducial parameters: for 
several values of a, models I and 2 have similar amplitudes 
of the MHD torque excess. Therefore, we expect the results 
presented in Section [4] to be representative of the migration 
of low mass planets when the interaction of the planet with 
the disc is linear. Additional complexity may arise when 
non-linear effects become important at larger planet masses 
and smaller softening lengths. 



lUribe et al.l (|201 lh also performed simulations of planet 
migration in 3D MHD turbulent discs. We do not attempt 
a detailed comparison with their results because the disc 
density profile and the strength of the turbulence signifi- 
cantly evolve during the time of their simulations. A qualita- 
tive comparison is however interesting. In their simulations 
of intermediate planet mass experiencing type I migration 
(with M p /(A/*/i 3 ) ~ 0.15 - 0.6), they observed a positive 
torque exerted on the planet, corresponding to outward mi- 
gration. They interpreted this result as a consequence of a 
strong and positive hydrodynamical corotation torque in a 
loc ally flat or increasin g surface density profile as described 
by iMasset et~aH (|2006r i. In light of our results, we suggest 
that the MHD torque excess might have played a prominent 
role in the positive torque they observed. Indeed, the den- 
sity and temperature profiles in their simulations (initially 
the same as Model 1, and later on showing a locally flat or 
increasing density profile) should lead to a positive MHD 
torque excess (see Figure [TT1 or equation I17[) . It is hard to 
make any prediction for the amplitude of the MHD torque 
excess given that they used different planet masses and disc 
aspect ratio than our fiducial values, but we note that the 
low value of the effective viscosity measured in their simu- 
lations (a ~ 1 - 2 x 10 -3 ) is favourable for a strong MHD 
torque excess. 



6 DISCUSSION AND CONCLUSION 
6.1 Summary 

Using two-dimensional numerical simulations of viscous re- 
sistive laminar discs, we have confirmed and explained the 
existence of an additional corotation torque due a weak az- 
imuthal magnetic field, which was observed in the 3D MHD 
turbulent disc simulations by BFNM11. It arises from the 
interaction between the horseshoe motion of the gas in the 
planet's coorbital region and an azimuthal magnetic field in 
the protoplanetary disc. We refer to this new torque com- 
ponent as the MHD torque excess. This is d i stinct from 
the magnetic resonances described bv lTerqu cm (2003), and 
arises in a different regime of magnetic field strength. In- 
deed, we considered magnetic fields that are weak enough 
not to prevent the horseshoe dynamics. This is the case in 
our simulations where the Alfven speed is smaller than the 
shear velocity at the separatrix of the horseshoe region. We 
do not see any evidence of the presence of magnetic reso- 
nances in this regime. As discussed in Section \2. 31 this may 
be due either to the horseshoe dynamics (the presence of the 
horseshoe trajectories radically changes the shear, such that 
the description in terms of magnetic resonances is not ap- 
propriate anymore) or the effect of resistive diffusion, which 
we included to describe the effect of turbulence. 

The physical origin of this new component of the coro- 
tation torque has been interpreted in the following way. The 
horseshoe motion of the gas near the planet tends to accu- 
mulate the magnetic flux contained in the horseshoe region 
on the downstream separatrix. This results in a region of 
enhanced magnetic pressure around the downstream separa- 
trices, where the gas surface density is reduced to maintain 
approximate total pressure balance. The MHD torque ex- 
cess most likely stem from an azimuthal asymmetry of this 
underdense region. This is linked with an azimuthal shift of 
the stagnation point with respect to the planet, since the un- 
derdensity is maximum close to the stagnation point. Such 
an asymmetry is always observed in the presence of a mag- 
netic field, and increases with the magnetic field strength. A 
positive azimuth of the stagnation point leads to a negative 
torque excess, while a negative azimuth of the stagnation 
point leads to a positive torque excess. This sign is governed 
by a hydrodynamic process determined by the temperature 
and density gradients. The amplitude of the torque excess 
on the other hand is governed by a MHD process determined 
by the magnetic field strength and diffusivity. 

The properties of the MHD torque excess may be sum- 
marised in the following way: 

• Its magnitude is roughly independent of the gradients 
of magnetic field strength, density and temperature, which 
have at most an effect of order 10 — 20%. In the non- 
dimensional problem, it therefore depends only on the mag- 
netic field strength (through /3), the diffusion coefficients 
(a and V), the planet-to-primary mass ratio M p /M*, the 
disc's aspect ratio h, and the ratio of softening length of 
the planet gravitational potential to the pressure scaleheight 
e/H. By varying the parameters a and /3, we obtained the 
following scaling of the amplitude of the MHD torque excess: 
1 oc l/(/3 2 a 4 ).Thusit increases very strongly when the mag- 
netic strength is increased or when the turbulent resistivity 
is decreased. The amplitude is also roughly consistent with 
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being inversely proportional to the softening length of the 
planet potential. The dependence on M p /M*, and h is left 
for future work. 

• The sign of the MHD torque excess is determined by 
whether the azimuth of the stagnation point is located ahead 
of or behind the planet, and this depends solely on the den- 
sity and temperature gradients at the planet's location. It is 
therefore independent of the strength and gradient of mag- 
netic field and of the viscosity and resistivity. Our numerical 
results show that the sign of the MHD torque is determined 
by the slope indices of the temperature profile q, and surface 
density profile p in the following manner: 

sign(r) = sign(p -1.3q + 0.6) (20) 

We compared the results of our 2D laminar simulations 
with the 3D turbulent simulations of BFNM11. We found 
good qualitative agreement for the sign of the MHD torque 
excess and for the torque density distribution, and their de- 
pendence on the density and temperature profiles. The quan- 
titative comparison of the MHD torque excess amplitude is 
complicated by its steep dependence on the effective resis- 
tivity. The turbulent simulations results can be reproduced 
with values of the magnetic Prandtl number ranging from 
0.5 to 2 depending on the disc model (density and temper- 
ature profiles). 

6.2 Discussion 

We have estimated the MHD torque excess by performing 
non linear numerical simulations and found that it can have 
a significant impact on the migration. The amplitude of the 
MHD torque excess was found to be a very steep function of 
the magnetic field strength and diffusivity (through a and 
/?), which are rather uncertain. The magnetic field strength 
in protoplanetary discs is so far poorly constrained observa- 
tionally, but the SPIRou spectropolarimeter might provide 
more information in the future. From a theoretical perspec- 
tive, numerical simulations of MHD turbulent discs suggest 
values of a between a few times 10 -3 and a fe w times 10~ 2 in 
the a bsence of a net vertical magnetic flux (|Fromang et al.l 
l2007h . This range of values is also compatible with the ob- 
served characteristic lifetime of protoplanetary discs of a 
few millions years. The fiducial parameters used in this pa- 
per (a = 5.10 -3 , ft = 100) are thus plausible for MRI active 
regions of protoplanetary discs, and they lead to a strong 
MHD torque excess that is able to reverse the direction 
of migration. This suggests that the migration of low mass 
planets in MRI active regions of protoplanetary discs could 
be significantly affected by the magnetic field. This is fur- 
ther suggested by the comparison with MHD turbulent sim- 
ulations of planet migration detailed in Section [SI although 
further studies are certainly necessary to better characterise 
the influence of the MHD torque excess in self-consistent 
simulations with MHD turbulence. 

Could the MHD torque excess also be significant in dead 
zones, where the ionisation is too low for MHD turbulence 
to be sustained ? This possibility may seem counterintu- 
itive at first sight because dead zones are characterised by 
a large physical resistivity that prevents turbulence, while 
a significant MHD torque excess needs the resistivity to be 
low enough. However both criteria may be satisfied at the 
same time. Indeed the critical magnetic Reynolds number 



below which the turbulence is shut off is thought to lie in 
the range 10 2 — 10 4 depen ding on the presence or not of a 
net vertical mag netic flux (|Fleming et alJl200Ch . which cor- 
responds to a/V = 10~ 4 — 10~ 2 . Our results show that the 
MHD torque excess can be significant for this range of resis- 
tivity if an azimuthal m agnetic field is presen t in the dead 
zone. The simulations of lTurner fc Sanol (|2008l ) suggest that 
such an azimuthal magnetic field can arise in a dead zone 
owing to the interaction with the active surface layers of the 
disc. The presence of a significant MHD torque excess tak- 
ing place in dead zones is therefore an intriguing possibility 
which should be further studied. 

What sign of the MHD torque excess should one expect 
for disc profiles typical of protoplanetary discs ? Answering 
this question is difficult because the structure of protoplan- 
etary discs is still poorly known especially in their inner 
parts. Due to angular resolution limitations, observations 
can so far only constrain the outer parts of protoplanetary 
discs at distances of several tens of AU. In this disc region, 
they suggest the following tem perature and density pro files: 
q e [-0.7,-0.4], p e [-1,0] | Williams fc Ciezall201ll . and 
references therein). This region of parameter space lies en- 
tirely in the region where the MHD torque excess is positive 
(see Figure [TTJ. We therefore expect that in outer regions of 
protoplanetary discs, the MHD torque excess can slow down 
or reverse an otherwise inwards type I migration. The den- 
sity and temperature profiles in the inner part of the disc are 
much more uncertain so far, though the ALMA interferom- 
eter should enable its study in the near future. One might 
speculate that very negative values of p are unlikely and 
therefore that the MHD torque excess is likely positive. We 
should caution that all this discussion is assuming smooth 
density and temperature profiles. However what really mat- 
ters is the local gradient near the planet in a region of radial 
width ~ H, which may be very different from the global 
ones in some special locations. For example, a locally flat or 
positive surface density profile has been argued to provide 
a planet trap thank s to a strong positive corotation torque 
l|Masset et alj 120061 ). We note that the MHD torque excess 
is positive for such a density profile (unless the temperature 
is strongly increasing at the same location) and could thus 
make the planet trap mechanism even more effective. 

As discussed above, an interesting consequence of the 
MHD torque excess is that it could slow down or re- 
verse an otherwise inward type I migration. This is es- 
pecially interesting in the outer parts of the disc be- 
cause they are radiatively efficie nt (typically at r > 15 AU 
IjPaardekooper fc MeUemall2008i )). so that the hydrodynam- 
ical corotation torque is probably small and thus inefficient 
at counteracting the inward migration. Reversing the migra- 
tion at large radii could be important to explain the observa- 
tions by direct imaging of planets at several te ns of AU from 
their star li ke the 4 planets system HR 8799 (|Marois et alj 
120081 . |2010| ). We speculate that the MHD torque excess 
could lead to an outward migration to large distances and/or 
might in some cases prevent the inward migration of plan- 
ets formed at large di stances, which were sh own to migrate 
inwards very fast by iBaruteau et alj (|20 1 lh . As shown in 
this article, the positive torque leading to outward migra- 
tion can be as large as the diff erential Lindblad tor que. Us- 
ing the standard expression of iTanaka et al.1 (|2002l ) for the 
differential Lindblad torque in 3D, and a minimum mass so- 
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lar nebula disc model, we find that a 10 Earth mass planet 
could migrate to 50 AU in 1.5 x 10 6 years. There is therefore 
enough time for the migration to occur during the lifetime of 
the protoplanetary disc. That the MHD torque excess could 
lead to a sustained outward migration over long distances 
remains, however, to be proven. Indeed, as the planet mi- 
grates outwards the ratio of horseshoe U-turn time to the 
diffusion time across the horseshoe region is likely to vary, 
which would change the efficiency of the MHD torque excess. 
In this respect, it would be important to investigate the de- 
pendence of the MHD torque excess on the disc aspect ratio 
and the planet mass. 

Finally, we stress that we have considered in this study 
only a simple magnetic configuration: a large scale purely az- 
imuthal magnetic field. In reality, the field geometry is likely 
to be much more complex. For turbulent regions of the disc, 
it should be checked further if the effects of a disordered field 
structure can be described by a mean field with diffusion co- 
efficients, as is assumed here. Another question is the effect 
of a vertical magnetic field. The strength of the vertical mag- 
netic field in accretion discs is still debated. It depends on 
the relative efficiency of the advection by the accretion flow 
and the diffusion by the tu rbulent resistivity, whic h is still an 
active subject of research l|Guilet fe Ogilvial2012l . and refer- 
ences therein). The presence of jets in many T-Tauri stars 
may suggest that a strong vertical magnetic field threads at 
least the inner part of protoplanetary discs according to the 
magneto- centrifugal model of jet launched from an ac cre- 
tion disc (|Blandford & Pavndfl982l ; iFerreira et alj|2006j ). It 
would thus also be important to investigate the influence of 
a vertical magnetic fi eld on planetary migration, as started 
bv lMuto et all (|200Sh . 
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APPENDIX : CONVERGENCE TEST AND 
CODE COMPARISON 

We have tested the convergence and compared the results 
with the two codes used in this article, RAMSES and NIR- 
VANA. For this purpose we chose to use Model 2 with the 
fiducial parameters both in MHD (/3 = 100) and in hy- 
drodynamics cases. We used three different resolutions: the 
nominal resolution used throughout this paper where the 
half width of the horseshoe region is resolved by about 9 
cells (NR, blue curves), a twice lower resolution (LR, black 
curves) and a twice higher resolution (HR, red curves). No- 
ticeable differences are observed at low resolution although 
the numerical error on the torque remains modest (5 — 10%): 
in particular the two peaks in the torque density distri- 
bution are slightly less pronounced than at higher resolu- 
tion. The agreement between the high and nominal resolu- 
tion is remarkable with the code RAMSES, as the torques 
and torque density distributions are almost indistinguish- 
able. This agreement is slightly less good for the code NIR- 



VANA but still quite acceptable with a 5% difference in the 
torque. Finally the results of the two codes agree reason- 
able well with each other, with a difference of ~ 5% in the 
torques computed at nominal or high resolutions. 
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Figure 20. Convergence test of the two codes RAMSES (left panels) and NIRVANA (right panels) for the model 2 (p = —3/2 and 
q = 0) with our fiducial parameters. The upper panels show the total torque exerted on the planet, MHD simulations being represented 
with full lines and hydrodynamical ones with dashed lines. The lower panels show the torque density distribution as a function of radius 
near the planet for the MHD simulations (normalised by the maximum value of the high resolution run). The vertical dashed line shows 
the corotation radius, while the dotted lines show the extent of the horseshoe region. Three resolutions have been used: the nominal 
resolution used in this paper (NR, blue curves: 512 X 1024 for RAMSES and 468 X 1248 for NIRVANA runs), a twice lower resolution 
(LR, black curves: 256 X 512 for RAMSES and 222 X 612 for NIRVANA runs), and a twice higher resolution (HR, red curves: 1024 X 2048 
for RAMSES and 936 X 2496 for NIRVANA runs). Nominal and high resolution arc in excellent agreement for the RAMSES runs, and 
reasonable agreement for the NIRVANA runs. The torques computed with both codes at nominal or high resolutions agree within ~ 5%. 
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